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Abstract 

For the radiative Bhabha calibration of BaBar's electromagnetic calorimeter, the 
measured energy of a photon cluster is being compared with the energy obtained via 
a kinematic fit involving other quantities from that event. The details of the fitting 
algorithm are described in this note, together with its derivation and checks that ensure 
that the fitting routine is working properly. 
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1 Introduction 



Radiative Bhabhas can be used as one of the calibrations of the BaBar electromagnetic calori- 
meter (EMC). Radiative Bhabha events (e~e + — > e~e + j) deposit photons over a large energy 
range everywhere in the calorimeter. If the momenta of the incoming and outgoing electrons and 
positrons, as well as the photon's angular position are known, the photon energy can be obtained 
via a kinematic fit. This fit results in an absolute measurement of the photon energy which then 
can be compared to the measured photon energy to obtain calibration constants. 

The radiative Bhabha module is part of BaBar's Online Prompt Reconstruction (OPR) exe- 
cutable. Initial cuts select good electrons, positrons, and photons. Then all possible combinations 
of triplets (one electron, one positron, one photon) are formed. Each triplet is sent to the fitting 
routine to calculate its Xest> the "estimated x 2 "- The triplet with the lowest xLt is then sub- 
mitted to the full kinematic fit which returns, among other quantities, the fitted photon energy 
Efa and the error matrix of the fitted quantities. The ratio E meSuS /Ef^ is later used to calibrate 
the calorimeter. Note that no information on the measured photon energy E^^s goes into the 
kinematic fit or xLt- 

This note is the complete documentation on the algorithm for fitting the radiative Bhabha 
events for the purpose of calibrating the calorimeter. It describes the whole fitting procedure: the 
quantities for the kinematic fit and Xesti the derivation and formulas for Xest'i the derivation and 
algorithm for the kinematic fit; tests to check the quality of the kinematic fit. The note details all 
formulas which go into the computer program so that the program can be checked directly against 
this document. The derivations contain more details than needed to understand the concept, but 
the details help to derive, check and recheck all necessary formulas. Actual results of the fitting 
procedure using real data are not included in this note to keep it a pure code documentation. 



2 Defining the quantities and constraints 
2.1 Measured quantities 

From the experiment come the following measurements, which shall form the 14-dimensional 



vector y: 
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The momenta of the incoming electron and positron and their errors are changing run-by-run. 
The errors of the incoming leptons are given as covariance matrices: 
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The errors on Pj_ and Pj + are assumed to be independent. 
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The errors on #„ 7 and 



'o 7 appear in the current analysis without ^-^-correlations since they were 
found to be negligibly small, but we still use this 2x2 sub-set of the larger 4x4 error matrix of 
the EmcCluster: 
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All the errors can be combined in one 14 x 14 error matrix V^n. Its format is like this: 
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2.2 Quantities for the kinematic fit 

The kinematic fit determines the following numbers: 
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The variables f\ to /14 have corresponding measurements. The variable hi, the photon energy, 
is called a "hidden variable". The vector a shall be defined as a 19-element composite of / 
(14 elements), h (1 element), and A (4 elements). 



2.3 Constraints 

We have four constraint equations that have to be satisfied in the kinematic fit: 

Pi x - + pi x+ — pox- — Pox+ — Ef 7 sin 9 cos </>/ 7 = momentum in x 

Piy- + Piy+ — POy- — P0y+ — Efj sin 9 j 7 sin 4>f 1 = momentum in y 

Piz- + Piz+ ~ Poz- - Poz+ - Efa cos 6*/ 7 = momentum in z 

Ei- + E i+ - E ^ - E 0+ - E }1 =0 energy 

Here we use, e.g., 
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y/f? + fi + fi + m\ 



3 The estimated \ 2 '- Xest 

This function is calculated for any given electron-positron-gamma triplet to determine which 
triplet should be used for the kinematic fit. At the end of this subsection, we will have a complete 
analytical formula for calculating Xest- 

The formula is based on the difference between the initial and final momentum, P = Pi — P Q . 
The initial momentum Pi is the sum of the momenta of the incoming electron and positron as 
defined earlier: Pj_ and Pj+. The measured momenta of the outgoing electron, positron are 
given by Po- an d Po+- 

For the outgoing photon, we only have its angles #o 7 and 4>oj. Using the energy constraint 

Ej = Ei^ + Ei + — i?o- — Eq + 
we may substitute the unknown photon energy E y with measured values, and we obtain: 
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Of course, n is the normal vector, the direction of the photon. 
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Calculating the difference to form vector P is easy: 
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In the ideal world, this vector would be exactly zero. For its error matrix V p , we convert V a \\, the 
error matrix of y, via a transformation matrix T into Vp: 
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For the transformation matrix T we have to calculate expressions like -M 
j = x,y,z: 

djEjUj) _ Pjx- 
q p. -p. n i 

ul ix— 1 i — 

The transformation matrix is a 3 x 14 matrix: 
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What is the meaning of this x 2 ? We can say that the 14 input variables are used to measure 
P, and Xest tells us the deviation of the measured P from the expected P, which is zero. 



4 The kinematic fit 

For the derivation of the kinematic fit algorithm, we follow the description of Louis Lyons, page 
151, 152 §. 

4.1 The x 2 -Function 

The real x 2 -function can be written down in the following way: 

X 2 = (f-myv^if-m) 

+ Ai \p xi - + p xi+ - p xo - - p xo+ - E 7 sin 7 cos </> 7 ] 

+ A 2 [Pyi- + Pyi+ - Pyo- ~ Pyo+ ~ Ej S1U 8 y Sill <j)y] 

+ A 3 [p zi - + p zi+ - p zo - - p zo+ - E y cos 6y] 
+ A4 [Ei- + Ei + — E a _ — E a+ — Ey] 

The constraint equations are here included via Lagrange multipliers. To minimize this x 2 -> 
we could use a standard package like MINUIT, but standard packages are always slower than 
specially adapted code. Since the ^-minimization is being done millions of times, it pays off to 
write special code for the minimization. In addition, MINUIT is not supported in BaBar's Online 
Prompt Reconstruction. 

4.2 Derivation of kinematic fit algorithm 

At the minimum of x 2 > its first derivatives are to be zero. Lyons uses for this the following 
equations: 



dx 2 

da.i 

dx 2 
dh 

dx 2 



The three equation sets can be written as:fj] 



for i = 1 to 14 
here h = E y = ais 
here Ai = aie etc. 



2G(f -m)+D t X = 
E t A = 
c = 



^The factor 2 in front of V^ 1 is missing in Lyons' book |]J. We could easily remove this factor from our formulas 
by re-defining the Lagrange multipliers in the y 2 -function with a factor 2. This would not change the fit result or 
errors, as long as the subsequent calculations were carried out consistently. 
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where G is the 14 x 14 inverse error matrix of the measurements which we also call V\\ ■ 
/dCi/dai ... dCi/dau\ 



D 



and 



E = 



dCijdax ... dC2/dai4 
dC^/dati ... dCs/dau 
XdC^/dai ... dC^/dau/ 

We now expand the constraint equations C around fo and ho, and we obtain for the four equations 
Ck with k = 1 to 4: 
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and we see: 

= Z 

This is the equation we have to solve. Since the constraint equations C = contain non-linear 
functions like sin# 7 , Eq. (|l]) is only an approximation, and we have to iterate as described in the 
next section. 



4.3 Recipe for the kinematic fit algorithm 

The matrix M and the vector Z are functions of the measurements and their error matrices as 
well as of the parameters ex. The vector Y is, as mentioned above, 

f — m 
h — ho 
A 



and can be calculated with: 



Y = M" 1 Z 



Here is the iteration: Initially, we will use for the fit quantities / = m, i.e. the measured 
quantities. For h = ho, we calculate the photon energy via simple energy conservation. These 
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together with the measured quantities allow us to calculate M and Z. We multiply the inverse 
of M with Z and obtain Y. This result will then give us a better set of / and h, which we again 
use to calculate M and Z, and then a better Y. And we continue until our constraint equations 
are sufficiently fulfilled and the quantities / and h are stable. 

It might happen that the iteration does not converge at the minimum, but wanders off into 
unphysical numbers. In that case, it would be good to have a certain boundary box around the 
point. If the step would make the point lie outside the box, then the program would change the 
step so that the point would be back inside. It might be good to implement this, although the 
radiative Bhabha fitting does not seem to need this part of the algorithm. 



4.4 Details of matrices and vectors used in the kinematic fit 

We define the following variables: 



Ei- 


for i = 


1) 2, 3 (j>i x —i Piy~j Viz—) 
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Ei 



g = j 1 for i < 6 (p ix -, Piy-, P iz -, P ix+ , Piy+, P iz+ ) 

\ -1 for i > 6 (jPQx—i P0y—7 POz—i PQx+j P0y+i P0z+) 

For 4 x 14 matrix D we need the following expressions: 
Row j = 1 to 3, columns i = 1 to 12: 

dCj _ J Si if i = j or i = j + 3 or i = j + 6 or i = j + 9 
dai 1 else 

Row j = 1, column i = 13: 

— = — -C//7 cos ttfy cos0j 7 = — a 15 cosai3 cosai4 

a«i3 

Row j = 1, column i = 14: 

dC\ 
dan 

Row j = 2, column i = 13: 

dC 2 
dai 3 

Row j = 2, column z = 14: 

dC 2 
dai4 

Row j = 3, column i = 13: 

dC 3 



= Ef^ sin^j 7 sin0j 7 = Q15 sin 013 sin 0:14 



-Ef^ cos^j 7 sin0j 7 = — a 15 cos«i3 sinai4 



Ef^ sin#j 7 cos0/ 7 = — «i5 sinai3 cosai4 



da 



13 



= Ef y sinOfy = Q15 sin«i3 
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Row j = 3, column i = 14: 
Row j = 4, columns i = 1 to 12: 
Row j = 4, column z = 13 and 14: 



= 



dC 3 
dai4 

dC^ ai 
den 1 Ei 



dd dC 4 



dai3 da 



= 
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The 4x1 matrix E is: 



E 



I — sin 9 / 7 cos (f)fj \ 
— sin 9 j 7 sin </> j 7 
— cos 9 j 7 
-1 



t 



( — sin ai3 cos ol\± \ 
— sin«i3 sinai4 

— COS «13 
-1 



4.5 The error matrix of the fit 

The second partial derivatives of \ 2 appear in the error matrix of the fit parameters: 



H 



i d 2 x 2 



2 dctidctj 



So in our case, H is a 19 x 19 matrix. The detailed expressions for the second derivatives of % 2 
will be given in the following section. 



4.6 Tests for goodness of fit 

After completing the iteration on the kinematic fit, one wants to make sure that all quantities are 
indeed correct. 

Besides the obvious tests that the constraint equations are satisfied, one can check that indeed 
a minimum was reached. For this, one may wiggle each final value a\ to «i4 and recalculate \ 2 - 
In our case we have in the x 2 -function the terms with the Lagrange multipliers. Just recalculating 
the x 2 function will not lead to correct results, since the found vector a is a minimum only when 
also requiring the constraints. So one has to redo the fit while forcing the selected element of a 
to the off-minimum value. 

This wiggling allows us to map out the minimum, and it also tells us whether the fit error 
returned for that parameter is reasonable. If we fix £j 7 to be ±lc 7 fi t way from the real fit result, 
then the \ 2 should rise by 1 in either direction. When mapping out this rise, one will see the shape 
of a parabola. When the formulas are complicated and/or one is far away from the minimum, the 
parabola will be distorted. 

In our case, we can indeed calculate the fit error for -E/ 7 , but if this would be impossible, 
one can find the fit error by mapping out the minimum with the above described re-fitting with 
fixed values. The ±l<r-error is then defined to be where x 2 is 1 um t above the minimum. As 
mentioned, this function may be distorted when far away from the minimum. A complicated 
X 2 -function might even distort the ±lc-area. In this case, one can take the minimum and two 



9 



points very close to it, fit a parabola through these three points, and take the sigma from that 
parabola as the error. 

The same process also works for the hidden parameter (fitted photon energy) , and we definitely 
have to re-fit since the fitted photon energy only appears in the constraints, where the Lagrange 
multipliers would influence the outcome. 

Here is how we have to modify the formulas for re-fitting: 

4.6.1 Re- fitting with fixed Ef 1 

We want to redo the fit with the photon energy fixed to Er x = Ef^ + e. To the x 2 -function, we 
add the term 

+ X (E 7 - -Efix) 2 

where X is a large number compared to the original \ 2 ■ If we now minimize this new x 2 -function, 
the additional term adds a large penalty to any deviation of E^ from E^ x . 

Going through the derivation again, we find the following places that have to be changed in 
the code: 

• First partial derivative ^L. f or z = 15 [for i = k] has the additional term U +2X(E^ — Efix)". 

• No change to second partial derivatives. 

• Matrix M has the additional term "+2X" at (15,15). This means that the (15,15)-element 
of M is no longer zero. 

• Vector Z has an additional term at position 15: 

( ° 

Z = -2X(h - E &x ) 

\ -R 

These are all necessary changes. The iteration should converge again, but this time always 
result in E 1 = E^ x for sufficiently large X. 

4.6.2 Re-fitting with fixed 

Let us now wiggle one of the measurement variables a\ to 014. When fixing to /& = /fcfi X) we 
add the term 

+ X(f k -f k&x ) 2 

to the x 2 -function. Again, X is a large number compared to the original \ 2 • The following changes 
have to be made in the formulas of the algorithm: 

• The first partial derivative d\ 2 /dati gets for i = k the additional term "+2X(f k — /fefbc)". 

• Again no change to second partial derivatives. 

• Matrix M gets at position (k,k) the additional term "+2X". 

• Vector Z has at position k the entry ll —2X(mk — fk&x)" • 
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4.6.3 Confidence Level 



If all errors of the measurements are nicely described by Gaussian distributions, and if all events 
are what we think they are, i.e., (in our case) radiative Bhabhas, then the y 2 values of the fits 
should be distributed like the ^-distribution for n = 3 (3 because out fit is a 3-constraint fit). 
Instead of looking at the y 2 distributions directly, it is easier to map the x 2 to a flat distribution 
with values between and 1. This value is then called the confidence level (C.L.) of the event. If 
the x 2 is really distributed as it should be, the confidence level will have a flat distribution. 
So we are looking for two things in the C.L. distribution: 

(1) Most of the region should have a flat distribution. If not, the errors used in the fit might be 
too large or too small. If the errors are underestimated, the y 2 will be larger than expected, and 
the confidence level distribution will be tilted downward (when going from to 1). Vice- versa, 
if the errors are overestimated, the C.L. distribution will be tilted upward. More information on 
the validity of errors might be obtained from the "pull" distributions described later. 

(2) A peak at zero indicates events that do not fulfill the kinematics of radiative Bhabhas 
at all. They will result in very large x 2 (=very small C.L., close to zero). These events can come 
from backgrounds or misidentified tracks. What can we do? We can improve our selection criteria. 
Or we can cut out all events belonging to that peak, taking only those events that are part of the 
flat distribution. A cut on the confidence level is, of course, equivalent to a cut on y 2 - 

4.6.4 The "Pull" 

For each measured variable, one can plot the so-called "pull" ^] or "normalized stretch val- 
ues" 1@: 

meas — fit 
pull p = —^^=^^= 

V°"mcas — Cfit 

The minus sign in the square root comes from the strong correlation between the measured and the 
fitted quantity, and "still puzzles many users" |||. If all measured errors were estimated correctly 
and the conditions for the fit were satisfied {e.g., the event was really a radiative Bhabha event), 
then the pull quantity will be distributed like a Gaussian centered at with a = 1. If an error is 
for example overestimated, the pull quantity will have a more narrow distribution. In this case, 
the confidence level should also be affected, displaying a tilt in its distribution. 

To check whether a systematic increase or decrease of one or more errors would improve the 
pull and/or the confidence level distributions, one can redo the whole analysis with increased or 
decreased errors. Perhaps one can find a set of corrections that create nice pull distributions and 
a nice confidence level distribution. If the errors are really not correct, one should talk with the 
colleagues who are responsible for the errors. However, abnormal pull quantities might not be 
always created by incorrect errors. Systematically shifted measurements could also cause such 
symptoms. 

5 x 2 -Function — First Derivatives 

For this set of equations, we will use the following notation: 

{Ai for i = 1, 4, 7, 10 (p ix -, p ix+ , po x -, Pox+) 

A 2 for i = 2, 5, 8, 11 {p iy -, p iy+ , p 0y ~, Po y +) 

A 3 for i = 3, 6, 9, 12 {p iz -, p iz+ , p 0z -, Poz+) 
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Now we calculate the first partial derivatives of the x 2 -function, i.e., the 19 equations dx 2 /dai 
For i = 1 to 12: 



3=1 



-q^.= 2 Y, Kuljtfj ~ m j) + s i L i + s ^Yi 



For i = 13: 



2 14 



= 2 ^ KiiHj(/j ~~ m i) ~~ ^1-^7 cos #7 cos 07 — A2-E 7 cos 6y sin 7 + X^E^ sin # 7 
O'Oj . =1 

14 

= 2 E F an'>i-™;) 

— ai6CKi5 cos «i3 cos «i4 — 017015 cos 013 sin ai4 + ai8Q;i5 sin 013 



For i = 14: 
<9x 



da. 



For i = 15: 



2 14 

= 2 ^ ^Kthlj (fj — m j) + X\E 1 sin Q 1 sin 7 — A2-E7 sin Q 1 cos 7 
j'=i 

14 

= 2 Kuij^^j - m j) + ai6"i5 sinQi 3 sinau - ai 7 ai 5 sinai 3 cosai 4 
3=1 



dx 

— — = — Ai sin^ 7 cos 7 — A2 sin 7 sin0 7 — A3 cos 6* 7 — A4 

= — «i6 sin «i3 cos ai4 — ai7 sin 013 sin a 14 — ais cos 013 — aig 



For i = 16: 



For i = 17: 



3a, ; 



For i = 18: 



dx 2 
da; 



p xi - + p xi+ - p x0 - - p x0 + - -E 7 sin Q 1 cos 7 
ai + 04 — 07 — aio — 015 sin ai3 cos ai4 

Pyi- + Py i+ - PyQ- ~ PyQ+ ~ ^ 7 Sin # 7 SO (j)^ 

a2 + a*, — a$ — an — 015 sin 013 sin ai4 

= Pzi- + Pzi+ ~ PzO- ~ PzO+ ~ E 1 cos 9 j 
= 03 + a$ - ag - 012 — ais cos ai3 



For i = 19: 



d x 2 

— — — + E i+ — £0- — -Eo+ — E 1 
oai 

= Ei- + Ei + — E'o- — -Eq+ — a i5 
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6 x 2- Fu nc tion — Second Derivatives 

For i = 1 to 12 and j = 1 to 12: 



2 Ku4'-^ A 4^ = 2^-^19 7 if = by definition 

all ij 



= 2V7,l else 



For i = 1 to 12 and j = 13 to 14: 

For i = 1 to 12 and j = 15: 

For i = 1 to 12 and j = 16 to 18: 

d 2 x 2 



d 2 X 2 i 



d 2 x 2 



ddjdai 



Si if Lj = Lj by definition 



dajdoLi % % 3 

= else 



For i = 1 to 12 and j = 19: 



d 2 X 2 U__ on 



dajdai Ei E; L 
For i = 13 and j = 13: 

— — - — = 2 V &lli - + Ai£^ 7 sin 6 y cos 7 + A2-E 7 sin 6* 7 sin 7 + X^Ej cos # 7 

(JfOLj(JfCxj r 

= 2 Kuij + a i6ai5 sin a± 3 cos a 14 + 017015 sin «i 3 sin a 14 + ai 8 a\ 5 cos a 13 
For i = 13 and j = 14: 

— — - — = 2V aUij ■ + Ai-E 7 cos # 7 sin</> 7 — A2-E' 7 cos 6> 7 cos </> 7 

C/Cxj (JijL 'i 

= 2 Kj]^- + CK16CK15 cos ai3 sin «i4 — 017015 cos 013 cos Q14 

For i = 13 and j = 15: 

2 

— — - — = — Ai cos 6* 7 cos 7 — A2 cos 7 sin0 7 + A3 sin6* 7 
oajOdi 

= — «i6 cos «i3 cos ai4 — 017 cos ai3 sin ai 4 + ot\% sin 013 



For i = 13 and j = 16: 

d 2 x 2 



= — E 1 COS Oy COS 7 = — «15 COS «13 COS «14 
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For i = 13 and j = 17: 

d 2 X 2 

q — q — = — cos # 7 sin 7 = — Q15 cos 013 sin ai4 

(JvjL j (JvjL^ 

For i = 13 and j = 18: 

d 2 X 2 

q — ^ — = i£ 7 sin 6> 7 = «i5 sin «i3 

CO, j kJLx^ 

For i = 13 and j = 19: 

— - — = o 

dajdai 

For i = 14 and j = 14: 

#V -1 
— — - — = 2 V alli ■ + A i E 1 sin # 7 cos </> 7 + A2 -E 7 sin 6> 7 sin 7 

CJCxj CJCXj / 

= 2 KiHj + a i6«i5 sin CK13 cos «i4 + 017015 sin ai 3 sin a 14 

For i = 14 and j = 15: 

«9V 

— — - — = Ai sin 6 y sin 7 — A2 sin # 7 cos </> 7 

CJCxj CJCX'l 

= a\6 sin «i3 sin ai4 — CK17 sin ce\3 cos «i4 

For i = 14 and j = 16: 

— — - — = i? 7 sin # 7 sin 7 

dajdai 

= ais sin«i3 sin«i4 



For i = 14 and j = 17: 

'.. - .".I 1 1 1/- I U> I .<- 



= --E , - v sin6U,cos< 



For i = 14 and j = 18 and 19: 
For i = 15 and j = 15: 
For i = 15 and j = 16: 

For i = 15 and j = 17: 



— «15 8111 «13 cos a 14 

d 2 x 2 



dajdai 







- 



dajdai 



dajdai 



sin # 7 cos <^> 7 = — sin 0:13 cos «i4 



dajdai 



— sin c/ 7 sm <p 7 = — sm Q13 sin ai4 
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For i = 15 and j = 18: 




= — cos 0~ = — COS «13 



For i = 15 and j = 19: 




-1 



For t = 16 to 19 and j = 16 to 19: 
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